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1  Introduction 

The  counterflow  configuration  provides  a  comprehensive  framework  for  studying  the  character¬ 
istics  of  non-premixed  laminar  and  turbulent  flame  problems  [1,  2,  3,  4,  5,  6].  However,  apart 
from  the  simplified  one-dimensional  spatial  models,  the  fidelity  of  direct  numerical  simulations 
(DNSs)  for  the  counterflow  configuration  in  terms  of  robustness/accuracy  exhibit  significant  sensi¬ 
tivity  to  the  boundary  condition  (BC)  treatment  [7,  8,  9] .  This  behavior  is  mainly  due  to  the  multi¬ 
directional  character  of  the  flow  at  the  boundaries  which  must  be  properly  accounted  by  the  BCs.  To 
mitigate  this  problem,  Yoo  et  al.  [7]  developed  improved  BCs  based  on  the  Navier-Stokes  Char¬ 
acteristic  Boundary  Conditions  [10]  (NSCBCs)  for  laminar  and  turbulent  counterflow  flames.  The 
major  improvement  consisted  in  introducing  the  (no-longer  negligible)  transverse  terms  into  the 
Locally  One  Dimensional  Inviscid  (LODI)  relations  in  order  to  capture  the  multi-dimensional  ef¬ 
fects  at  the  inflow/outflow  boundaries.  However,  two  major  shortcomings  can  also  be  identified 
with  the  improved  NSCBCs.  First,  they  preserve  the  use  of  relaxation  coefficients  entering  the 
improved  LODI  relations  and  these  coefficients  are  problem  specific.  These  coefficients  provide 
an  optimal  balance  between  maintaining  the  prescribed  upstream  values  for  the  inflow  variables 
and  reducing  spurious  wave  reflections,  and  usually  must  be  determined  through  a  trial  and  error 
process;  this  is  an  expensive,  time-consuming  procedure.  Second,  once  derived  for  real-gas,  the 
implementation  of  these  revised  NSCBCs  adds  a  considerable  computational  cost. 

In  this  work,  we  adopt  a  totally  different  approach  for  constructing  BCs  for  the  counterflow 
configuration  that  results  in  a  very  concise  framework.  In  particular,  we  combine  the  steady-state 
Navier-Stokes  equations  with  an  initial  flow  of  potential  type  to  construct  numerically  consistent 


1 


Sub  Topic:  Laminar  Flames 


algebraic  BCs  at  the  inflow  boundaries.  The  structure  of  the  paper  is  as  follows.  First,  the  gov¬ 
erning  equations  of  the  problem  are  presented.  A  detailed  discussion  on  the  construction  of  the 
initial  profile  of  the  flow  follows.  Then,  the  new  inflow  and  outflow  boundary  conditions  are  pro¬ 
vided.  The  specifics  regarding  the  numerical  implementation  follows.  Next,  the  validation  of  the 
proposed  numerical  differential-algebraic  BCs  is  performed.  Finally,  a  summary  is  given. 


2  Governing  Equations 

2,1  The  Navier-Stokes  equations 

The  Navier-Stokes  (NS)  equations  for  a  compressible  reacting  multicomponent  mixture  of  N 
species  expressed  in  terms  of  the  conservative  variables  U  =  [Um\  (m  =  1, ...  ,N  +  4)  can  be 
cast  in  compact  form  as 


d\J  dFl 
dt  dxi 


(1) 


where  F'  =  [F* J  (m  =  1, . . . ,  TV  +  4)  denotes  the  flux  vector  of  the  conservative  variables  along 
the  i- th  direction,  i.e. 


F*  =  \piii,  puxui  +  p8 u,  pu2Ui  +  p52i ,  pu3Ui  +  pS3i,  pipYu  . . . ,  pujYN_u  (pet  +  p)ui]r,  (2) 


and  the  vector  C  =  [C7m]  (m  =  1, . . . ,  N  +  4)  is  defined  as 


c  ^  [0.  9t«  dT2i 


dr 


3  j 


dJ 


i  j 


dxj  ’ 


dxj  ’  dxj  ’ 


dxj 


UJi, 


dJN-ij  .  dqj  dijijUi)  T 

- 7^ - r  U1n—1j  - 1  7. - 

OXj  OXj  OXj 


(3) 


The  usual  notation  applies  in  that  the  indices  i  and  j  follow  the  summation  convention,  t  is  time, 
Xj  is  the  /-th  spatial  coordinate,  p  is  the  mass  density  of  the  fluid,  Uj  is  the  j- th  component  of 
the  velocity  vector,  Yn  is  the  mass  fraction  of  species  n,  p  is  the  pressure,  ef  =  e  +  ^upp  is  the 
total  specific  energy  (e  denotes  the  specific  internal  energy),  Jnj  is  the  j- th  component  of  the  mass 
flux  vector  .Jr,  of  species  n,  qj  is  the  j -th  component  of  the  heat  flux  vector  q  and  un  is  the  mass 
reaction  rate  of  species  n.  Finally,  rl3  is  the  Newtonian  viscous  stress  tensor 


Tij  =  /i 


dui 

dxi 


2  duk  \ 

3  dxk  ij)  ’ 


(4) 


where  p  is  the  dynamic  viscosity  of  the  mixture  and  is  the  Kronecker  delta.  Body  forces  have 
been  neglected  in  the  present  equations. 

The  species-mass  and  heat  fluxes  have  been  given  in  detail  elsewhere  [11]  and  their  expressions 
are  based  on  the  fluctuation-dissipation  theory  [12,  13,  14].  The  mixing  rules  employed  for  the 
computation  of  the  mass  diffusion  coefficients  and  thermal  diffusion  factors  have  been  derived  by 
Harstad  and  Bellan  [15]  and  are  based  on  the  conjunction  of  nonequilibrium  thermodynamics  [12] 
in  the  limit  of  low  pressure  and  Grad’s  13-moment  theory  [16]. 

2,2  Transport  properties 

Transport  properties  were  computed  according  to  the  most  up-to-date  methods  [15,  17,  18,  19],  an 
example  being  provided  in  Masi  et  al.  [11]. 
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Figure  1 :  The  initial  potential  flow  entering  through  both  the  left  AB  and  right  A'B'  inlet  bound¬ 
aries,  and  exiting  through  the  AA'  and  BB'  boundaries.  The  solid  lines  represent  the  streamlines 
and  the  arrows  show  the  direction  of  the  flow. 


2.3  Equation  of  state 

All  NS  equations  are  coupled  with  the  Peng-Robinson  (PR)  equation  of  state  (EOS)  [20] 

RuT  &mix 

p  = - i - 2 — (5) 

from  which  T  and  p  are  obtained  as  an  iterative  solution  of  two  nonlinear  equations  that  satisfy 
both  the  values  of  p  and  e  as  obtained  from  the  solution  of  the  conservation  equations.  In  the  PR 
EOS,  v ph  is  the  PR  molar  volume  and  the  molar  volume  v  =  vpp  +  vs  where  vs  is  the  volume 
shift  computed  to  improve  the  accuracy  of  the  PR  EOS  for  high-pressure  conditions  [21],  while 
the  terms  amix  and  brmx  are  functions  of  T  and  X,  [11]. 


3  Initial  profile:  Inviscid  incompressible  potential  flow 

The  initial  profiles  correspond  to  a  two-dimensional  potential  flow  in  which  two  opposing  streams, 
each  of  constant  density  and  composition,  are  injected  from  the  two  boundaries  AB  and  A'B' 
of  length  Ly  in  the  ./  -direction;  the  streams  exit  the  domain  through  the  two  boundaries  AB  and 
A'B'  of  length  Lx  at  the  //-direction  as  illustrated  in  Fig.  1 .  Reference  quantities  (subscript  ref ) 
constant  in  time  are  first  specified  at  the  left  (superscript  /)  and  the  right  (superscript  r)  boundary 
points  C  and  C' .  The  reference  pressures  at  the  two  inlets  are  set  plref  =  plref  =  pref,  so  that  for 
the  initial  flow,  the  stagnation  point  (S)  is  always  located  at  the  centerline  with  x$  =  Lx/ 2  and  also 
ys 

The  initial  density  and  mass  fraction  fields  are  given  as 


r  0/  \  v^O /  M  f  [plrefiYnref\  if  0<X<Lx/2,Wy 

b(*,v),n(*.v)]  =  (fr£>y£]  if  Lx/2  <  x  <  Lx,Vy  ’ 


(6) 


where  a  discontinuity  arises  in  p°(x:y)  and  Yf(x1  y)  at  x  =  Lx/2,\/y  when  plrf,j  prr(,j  and 
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Y'i  ref  7^  y^ref  •  The  components  of  the  initial  velocity  field  are 


,oc  _  2k(o;  Lx/2)  ,  rw  „o  _  “2k{x  Lx/ 2) 


«“(*>?/)  = 


sjplef/p', 


0<x  <  Lx/ 2;  w°(x,2/)  = 


r 

re/ 


^ Pref/  Pref 


>  -fix/2  x  Y  Lx 


,0f„  _  ‘^,K{y  Ly/ 2)  ,  n/  Ly/ 2) 


w«(®,?/)  = 


r 
ref 


0  <  x  <  Lx/ 2;  uj(®,2/)  = 


1  +  \j  Pref /Plref 


(?) 

j  -Tx/2  <  x  <  Lx  (8) 


where  /c  is  the  strain  rate.  ux(x,y)  is  continuous  \/(x,y),  while  uy(x,y)  has  a  discontinuity  at 

x  =  Lx/ 2,  \/y  when  pj.e/  7^  p^e/. 

Using  (i)  the  common  reference  pressure,  pre/,  of  points  G'  and  C' ,  (ii)  the  Bernoulli  equation 
which  holds  for  inviscid  incompressible  flows  and  (iii)  the  fact  that  for  irrotational  (i.e.  poten¬ 
tial)  flows  the  total  pressure,  as  given  by  the  Bernoulli  equation,  is  the  same  for  all  points  of  the 
flow  (otherwise  the  Bernoulli  equation  would  only  hold  along  a  streamline),  we  obtain  the  initial 
pressure  field,  p°,  as 


P°(x,y)  =  Pref  +  -plref  ( U°x(xC,yC )2  +  Vc))‘ 


r\  Pref 


(u°x(x,y))2  +  (iiy(x,  y)Y 


for  0  <  x  <  -fix/2,  \/y,  where  xc  =  0  and  yc  =  Ly/2,  while  for  Lx/2  <  x  <  Lx ,  My 


9) 


p°(x,y)  =PrefJr~Pref  {Ul  {^C' ,  VC' ))  2  +  (®C' ,  VC')) ‘ 


1  r 

9  Pref 


K(^,i/))2  +  K  (x,y))' 


(10) 


where  Xc  =  Lx  and  yc  =  -fiy/2. 

By  construction,  the  initial  potential  flow  satisfies  the  steady-state  incompressible  Euler’s  equa¬ 
tions.  The  availability  of  p°(x,  y),p°(x ,  y),  Y®(x,  y), . . . ,  Y$(x,  y)  in  conjunction  with  the  EOS 
allows  the  determination  of  the  initial  temperature  field  of  the  flow,  T{)(x.  y).  With  the  calculated 
T°(x,y),  the  initial  specific  internal  energy  field  is  obtained.  It  follows  that  the  initial  field  of 
conservatives  variables,  U°(x,  y ),  can  then  be  readily  constructed. 

Finally,  we  note  that  for  the  counterflow  configuration  the  use  of  a  potential  flow  as  initial 
flow  is  consistent  with  the  fact  that  even  when  starting  computations  with  a  plug  flow  (i.e.  zero 
^/-components  of  the  velocity  vector  at  the  inlets),  the  flow  quickly  converges  to  a  potential  type. 


4  Boundary  conditions 

4.1  Inflow  boundary  conditions 

Instead  of  explicitly  imposing  timewise  constant  boundary  values  (also  known  as  hard  BCs)  for  the 
conservative  variables  at  the  two  inlets  AB  and  A'B',  we  impose  the  steady-state  NS  compressible 
form  of  the  conservation  equations  for  p,  pux ,  puy,  p\'\ , . . . ,  pYJv-i-  Using  these  equations  as  BCs 
means  that  the  initial  values  imposed  on  the  inlets  AB  and  A'B'  for  these  variables  are  implicitly 
constrained  to  be  constant  with  time.  Most  importantly,  the  transverse  terms  that  give  rise  to  multi¬ 
directional  effects  at  the  inflow  boundaries  are  inherently  present  in  the  BCs.  In  addition  to  the 
above  BCs,  the  thermodynamic  computation  of  the  energy  pet  =  pet(T(x,y),p(x,y),Yi(x,y), , 
Yn(x,  y))  is  used  for  the  boundary  values  of  pet  at  the  two  inlets,  where  et  is  calculated  using  the 
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current  values  of  p,  pux,  puy,  pY\ , . . . ,  pY^~  1  and  T  is  obtained  from  the  EOS,  so  that  the  values 
of  pet  are  readily  available. 

By  construction,  the  initial  values  of  p,  pux,  puy,  pY\ , . . . ,  pYN_]  at  the  inlets  (see  Section  3) 
satisfy  the  Euler  instead  of  the  NS  steady-state  equations.  Since  the  inlet  boundaries  are  usually 
located  away  from  the  flame,  these  initial  values  provide  a  good  initial  guess  for  the  NS  steady-state 
equations  which  can  be  corrected  [22]  to  account  the  viscous  and  species  mass-flux  effects. 

4,2  Outflow  boundary  conditions 

At  the  outflow  boundaries  A'  A  and  BB' ,  the  solution  satisfies  the  unsteady  compressible,  reactive 
flow  NS  equations  since  the  conditions  there  are  the  outcome  of  the  processes  taking  place  in  the 
interior  of  the  domain.  At  these  boundaries  we  assume  that  p  stays  constant  with  time,  taking  the 
value  initially  provided  by  the  Bernoulli  equation. 

5  Numerical  scheme 

An  eight-order  explicit  finite-difference  scheme  is  used  for  the  spatial  derivatives.  After  spatial 
discretization,  the  steady  NS  equations  become  algebraic  equations  at  a  node,  whereas  the  unsteady 
NS  equations  become  ordinary  differential  equations  at  a  node.  As  a  result,  at  the  nodes  of  the 
inflow  boundaries  AB  and  A'B',  the  BCs  are  numerically  of  algebraic  type.  In  addition,  auxiliary 
ghost  points  can  be  used  from  the  left  of  the  boundary  AB  and  the  right  of  the  boundary  A'B' 
(since  the  solution  is  known  there)  to  enhance  the  finite  difference  approximation  of  the  spatial 
derivatives  in  those  equations.  At  the  nodes  of  outflow  boundaries  A1  A  and  BB1,  the  BCs  are  of 
differential  type  and  one-sided  finite  differencing  is  used. 

The  resulting  numerical  system  of  differential-algebraic  equations  [23]  is  integrated  in  time 
with  the  differential-algebraic  solver  IDA  of  the  SUNDIALS  suite  [24] .  The  integration  method 
used  in  IDA  is  a  variable-order,  variable-coefficient  BDF  (Backward  Differentiation  Formula),  in 
fixed-leading-coefBcient  form  where  the  order  of  the  method  varies  between  1  and  5.  The  BDF 
method  can  handle  the  stiffness  introduced  in  the  numerical  integration  of  the  NS  equations  due  to 
the  presence  of  chemical  source  terms.  IDA  supports  parallel  computations  through  the  message 
passing  interface  (MP1)  protocol  and  provides  a  routine  that  computes  consistent  initial  conditions 
from  a  users’  initial  guess  [22,  25]. 

6  Reproducing  potential  flows 

To  validate  the  proposed  differential-algebraic  BCs  of  Section  4,  we  consider  two  symmetric  non¬ 
reacting  potential  flows  where  either  H2  or  02  enters  from  both  inlets.  The  goal  is  to  numer¬ 
ically  reproduce  these  two  steady-state  potential  flows  at  computational  times  characteristic  to 
reaction/diffusion  problems. 

Table  1:  Inflow  conditions  at  C  and  C1  (see  the  right  part  of  Fig.  1)  for  two  H2  and  one  02 
simulations.  Nx  and  Ny  are  the  number  of  points  in  the  x  and  y  directions.  Lx=  Ly=  10  mm. 


Case 

Species 

Pref  (kg/m3) 

Tref  (K) 

Pref  (bar) 

K  (S  *) 

Nx 

Ny 

A 

h2 

0.08077 

300 

1 

2  x  103 

104 

104 

B 

h2 

0.08077 

300 

1 

4  x  103 

144 

144 

- 

o2 

1.2837 

300 

1 

2  x  103 

200 

200 
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t  [ms]  t  [ms] 


Figure  2:  Temporal  variation  of  p.  vx,  vy,  p  and  YH.,  at  inlet  boundary  points  for  Case  A  of  Table 
1.  The  labeling  of  the  points  refers  to  the  right  part  of  Fig.  1. 

6.1  The  H2  symmetric  potential  flow 

For  the  symmetric  FI2  potential  flow,  the  two  situations  considered  differ  in  the  imposed  strain 
rates  (see  Table  1).  Due  to  space  constraints,  we  only  show  below  the  temporal  variation  of  the 
inlet  boundary  points  for  case  A  (Fig.  2)  and  the  spatial  variation  of  primitive  variables  at  t  =  6 
ms  for  case  B  (Fig.  3).  Figure  2  shows  the  temporal  evolution  of  the  boundary  values  of  p,  vx,  vy 
and  Yh2  at  points  on  the  left  inlet  boundary;  the  prescribed  initial  boundary  values  are  excellently 
maintained  by  the  imposed  NS  steady-state  BCs  at  the  two  inlets  AB  and  A'B' .  Figure  3  illustrates 
for  case  B  the  initial  (t  =  0  ms)  spatial  variation  of  p,vx,vy,p  and  YH.2  for  ^-sections  of  the 
computational  domain  and  the  spatial  variation  of  the  same  variables  obtained  at  t  =  6  ms.  Clearly, 
the  steady-state  potential  flow  is  accurately  reproduced  numerically  in  all  respects  for  the  larger  k 
value. 

6.2  The  02  symmetric  potential  flow 

For  the  symmetric  02  potential  flow  the  inflow  conditions  and  computational  mesh  are  given  on 
Table  1 .  The  results  of  the  simulation  for  02  represent  a  simulation  performed  with  a  fluid  having 
a  density  of  0(10)  larger  compared  to  that  of  H2  used  in  the  previous  two  simulations.  The  results 
are  not  illustrated,  but  they  show  the  same  high  fidelity  in  maintaining  the  imposed  initial  potential 
flow  as  in  the  H2  simulations. 

7  Summary  and  conclusions 

New  BCs  have  been  developed  for  the  counterflow  configuration  that  provide  a  very  concise  frame¬ 
work  inherently  accounting  for  the  multi-directional  character  of  the  flow  at  the  boundaries.  Upon 
discretization  of  the  steady-state  NS  equations,  the  BCs  at  the  inflow  boundaries  are  of  algebraic 
type,  whereas  upon  discretization  of  the  unsteady  NS  equations,  those  at  the  outflow  boundaries 
are  of  differential  type.  This  formulation  of  the  numerical  BCs  as  differential-algebraic  ones  re- 
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quires  the  use  of  a  numerical  integration  software  capable  of  handling  initial-value  problems  for 
differential-algebraic  systems  of  equations  [24] . 

For  the  non-reactive-flow  validation  of  the  differential-algebraic  BCs,  two  symmetric  non¬ 
reacting  potential  flows  were  considered  with  either  H2  or  02  as  injected  species  from  both  in¬ 
lets.  For  H2,  two  different  strain  rates  cases  where  examined  whereas  for  02  only  one  strain  rate 
flow  is  simulated.  Irrespective  of  the  species  or  the  strain  rate,  the  initial  potential  flow  was  re¬ 
produced  with  high  fidelity  over  the  entire  computational  domain  without  numerical  artifacts,  i.e. 
without  use  of  dissipative  filters  or  introduction  of  relaxation  constants  into  the  BCs. 
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